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The most massive black holes in our Universe form binaries at the 

centre of merging galaxies. The recent evidence for a gravitational- 

wave (GW) background from pulsar timing may constitute the first 
observation that these supermassive black-hole binaries (SMBHBs) 
merge. Yet, the most massive SMBHBs are out of reach of interferometric 
GW detectors and are exceedingly difficult to resolve individually 

with pulsar timing. These limitations call for unexplored strategies to 
detect individual SMBHBs in the uncharted frequency band <10° Hz to 
establish their abundance and decipher the coevolution with 

their host galaxies. Here we show that SMBHBs imprint detectable 
long-term modulations on GWs from stellar-mass binaries residing in 

the same galaxy at a distanced < 1 kpc. We determine that proposed 
decihertz GW interferometers sensitive to numerous stellar-mass binaries 
could uncover modulations from ~O(10 '-10*) SMBHBs with masses 
~O(10’-10*) M, out to redshift z ~ 3.5. This offers a unique opportunity to 
map the population of SMBHBs through cosmic time, which might remain 
inaccessible otherwise. 


® Check for updates 


We consider the situation in which we directly detect a GW signal a monochromatic modulation that can be written as a sinusoid 
from a stellar-mass compact binary at luminosity distance Dthatis (cf. equation (14)): 

accompanied by an SMBHB in its proximity at d < D, for example, 

at the centre of its galaxy (Fig. 1). The GWs emitted by the SMBHB CCE) = Amoa COS(2TY moat + Pmod)s (2) 


perturb spacetime at the location of the compact binary, which where fmoais the frequency of the modulating GW and its amplitude 
induces a small frequency modulation in the detected signal 


(cf. equation (11)): Amoa is of the order ~ (fmo) MSB d , With M moa being the chirp 

mass of the SMBHB. Thus, the detector receives modulated GWs from 

Ímeasurea (À © KO -fOO + OP), (D) the compact binary (cf. equation (21)), where the modulation is fully 

described by three additional physical parameters that are determined 

where f(t) is the unmodulated frequency of the compact binaryand by the SMBHB: an overall amplitude Amoa, a modulating frequency 
Z(t) <1is a small modulation. At the lowest order, the SMBHB causes fmoaand a phase @moa- 
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Fig. 1| Cartoon of the method proposed in this work. The presence of an 
SMBHB emitting GWs causes frequency modulations in the GW emission of a 
compact binary source at a distance d. The modulations can be observed over 
along observation time 7 with proposed decihertz GW detectors, at a distance 
D > d. We show how this scenario would allow decihertz detectors to indirectly 
probe the existence of SMBHBs within the -10’-10° M, mass range. 


To investigate how well these additional parameters can be meas- 
ured, we employ two methods. First, we calculate the Fisher matrix of 
the modulated post-Newtonian (PN) waveform of the compact binary 
that includes both higher-order modes and effective spin (Methods). 
The Fisher matrix approach allows us to efficiently survey a large part 
of the parameter space and to estimate the expected error of fitting 
the parameters of the sinusoidal modulation to the noisy data of a 
detector (Methods). Second, we run a series of full Bayesian parameter 
recovery tests employing numerical Markov chain Monte Carlo 
(MCMC) methods at a representative redshift of z = 0.84. In the latter 
method, we explore a large grid of injected values of A,,,gand finog tO 
find the curve delimiting detectable and undetectable modulations, 
where we consider a modulation to be detectable ifthe one-sigma width 
of the posterior distribution of Amoais smaller than the injected value. 
We find that this ensures that the detectable modulations are distin- 
guishable from the null hypothesis of no modulation (Amoa = 0) being 
present with extremely high confidence (see Methods for a more thor- 
ough discussion). Finally, as seen in Fig. 2, we can match the curve 
defined by this strict criterion to a corresponding relative error 


MAmoal!Amod = V (AAmod))/Amoa $107"? that results from the 
simpler Fisher matrix analysis. For the purposes of numerical effi- 


ciency, we use the scaling with redshift of the latter to extrapolate the 
stricter MCMC criterion. 

Modulations caused by SMBHBs can be most easily identified with 
a detector that can observe a large number of compact binaries witha 
high signal-to-noise ratio (SNR) and over a long observation time 7, 
which is limited by the mission lifetime of the detector. We find that 
these constraints single out decihertz GW detectors’ and the large 
populations of binary black holes (BBHs) and binary neutron stars 
(BNSs) that they are expected to observe. Thus, we adopt the sensitivity 
curves of the proposed space-based laser-interferometric GW detec- 
tors DECIGO and BBO”? as examples for any sufficiently sensitive detec- 
tor in the decihertz regime. In Fig. 2, we showcase the results of our 
Fisher matrix and MCMC analysis. As an example, we assume atypical 
BNS and BBH located at a representative redshift z = 0.84 and use the 
design sensitivity of DECIGO. The left panel shows the relative uncer- 
tainty AA moa/Amoa as a function of the modulating SMBHB frequency 
fmoa: The sensitivity curves have a characteristic shape that is similar 
to those of pulsar timing arrays*”, with a peak sensitivity around 
Soa © 1/T. The drop in sensitivity below finog S 1/T reflects the fact 
that the observation time needs to cover at least one period of the 
SMBHB GW iin order to establish its existence‘. Conversely, the sensitiv- 
ity degrades for higher frequencies following the fmoaT > 1 
limit of the hypergeometric functions in equation (19). The right 
panel of Fig. 2 shows the minimum detectable chirp mass 
Mmod = (Amod) Tfnoa) of the modulating SMBHB. Masses as 
low aS Mmoa ~ O(107) Mo can be detected if the compact binaries are 
located at a distance d = 1 pc to the SMBHB, as suggested by several 
binary formation channels (see below). 
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Fig. 2 | Detectability contours of the modulating SMBHB amplitude Amod 
experienced by a stellar-mass compact binary. The binaries are located in a 
fiducial host galaxy at redshift z = 0.84. The modulated GW signal from the 
compact binary is observed with DECIGO over T= 10 years. Left: to compute this 
sensitivity we vary the modulating amplitude Amoa and frequency fmoaofthe 
SMBHBand indicate the relative uncertainty AA mod/ Amoa = 10° and 10%, 
respectively, by which a given amplitude could be recovered with dashed contour 
lines, using parameter estimation with the Fisher matrix. Black contour lines 
assume that a GW150914 (ref. 100)-like BBH with a chirp mass of M = 28.0 Mais 
observed. Blue contour lines assume a GW170817 (ref. 68)-like BNS 

(M = 1.188 Mo). We also show for the BNS the sensitivity curve from a full 
Bayesian analysis with an MCMC (solid blue line) that corresponds to the Fisher 
matrix uncertainty AAmoaq/Amoa ~ 107"? (Methods and Supplementary Fig. 2). 
Right: we show the sensitivity in terms of the minimum detectable chirp mass 

M mod = (Amod) (ttf mod ys of the SMBHB for a fiducial distance d =1 pc to the 
BBH/BNS. The red line indicates the frequency above which our assumption of a 
monochromatic SMBHB breaks down. Note that the number of SMBHBs near this 
limit is anyways highly suppressed (Methods). 


In Fig. 3, we compare the sensitivity of our method with other 
currently operating and planned GW detectors. Since the SMBHBs 
perturb compact binaries in their proximity, we consider an equiva- 
lent strain sensitivity that corresponds to the strain amplitude the 
distant SMBHBs would produce on Earth. The resulting equivalent 
sensitivities of DECIGO/BBO in the nanohertz band are comparable to 
or better than those of current and planned pulsar timing arrays. For 
distances between the stellar-mass compact binary and the SMBHB 
of d=1pc, suggested by various compact object binary formation 
channels (see below), the obtained sensitivity would outperform the 
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Fig. 3 | Landscape of different GW detectors compared to the method 
proposed in this work. For each detector, we plot the dimensionless strain 
amplitude A = ViSn (solid lines) as function of frequency f (refs. 2,3,101-106), 
where the asterisks indicate currently operating detectors (in contrast to 
planned detectors). The black and blue lines show the equivalent strain 
sensitivity of our method, which corresponds to the GW strain amplitude the 
SMBHB would produce on Earth if it were indirectly detected (with 

AA mod/ Amoa = 10-2) through the modulations of a directly detected signal 
froma typical BBH (for example, GW150914 (ref. 100)) anda typical BNS (for 


10° 

f (Hz) 
example, GW170817 (ref. 68)), respectively. The dashed lines correspond to an 
observation with DECIGO and the dotted lines to an observation with BBO. We 
emphasize that this method can only detect SMBHBs that are concurrent witha 
compact object merger in its proximity at a given separation d, which limits the 
number of potential sources (see main text). Nevertheless, the resulting 
sensitivity can greatly outperform other nanohertz observatories, such as the 
European Pulsar Timing Array (EPTA), the International Pulsar Timing Array 
(IPTA) and SKA. 


anticipated sensitivity of the Square Kilometre Array (SKA) by two 
orders of magnitude. However, we stress that this striking sensitivity 
only applies to SMBHBs that are accompanied by a compact binary 
in their proximity, inherently limiting the number of detectable sys- 
tems. The sensitivity differences between BNSs and BBHs are caused 
by their different masses and frequency evolution (Methods). They 
are negligible for the purpose of estimating the detection rate (see 
below). Thus, we will henceforth refer to both as generic compact 
binary sources. 

Given the sensitivity curves above, we estimate the expected num- 
ber and properties of SMBHBs that cause detectable modulations on 
the GWs from compact binaries observed with decihertz instruments. 
For this purpose, we adopt an SMBHB population model based on the 
Millennium Simulation’. We then distribute compact binaries within 
galaxies following the latter’s total stellar mass and evaluate the prob- 
ability that they coincide with an SMBHB that produces a detectable 
modulation (Methods). 

Figure 4 shows that the distribution of detectable SMBHBs 
strongly peaks within the range of 10’-10° M,, which represents the 
best compromise between the total quantity of available SMBHBs 
and the strength of the modulation they produce. Most detections are 
limited to relatively low redshifts of 0.5 < z < 1. However, a substantial 
fraction of the potential detections is distributed in the large cosmo- 
logical volume enclosed between z= 1 and z = 3.5. We highlight that 
for the majority of detected compact binaries the expected distance 
measurement and angular resolution of DECIGO/BBO would be suf- 
ficient to identify the host galaxy of the SMBHB*, allowing for targeted 
multi-messenger searches of subparsec SMBHBs. This shows how our 
method can unlock individual GW-based detections of the most mas- 
sive SMBHBs in our Universe complementing current pulsar timing 
arrays and the upcoming LISA detector’ at higher redshifts and wider 
SMBHB separations, respectively. 

Figure 4 also shows the expected total number of detectable SMB- 
HBs as a function of the distance d and the compact binary merger 
rate R. The distance d at which the compact binaries typically reside 
fromthe centre of their galaxy is a consequence of how it was formed; 
multiple formation channels have been proposed and are being actively 
discussed in literature to explain the origin of the compact binary 
mergers observed with LIGO-Virgo-Kagra’’. For instance, in the field 


scenario compact binaries are formed far away from their galactic 
centre from the evolution of isolated massive stellar binaries, triples 
or more complex stellar systems”. In contrast, binary mergers near 
the galactic centre could be promoted by strong environmental effects 
in, for example, nuclear star clusters (NSCs)? or inside the disks of 
active galactic nuclei (AGNs)°*° °°. While the latter have only been inves- 
tigated for galactic centres that harbour a single supermassive black 
hole (SMBH), it is plausible to assume that they work similarly in galax- 
ies that host SMBHBs. Given the current LIGO-Virgo-Kagra observa- 
tions, it is uncertain which of these or other formation channels (for 
example, in globular clusters or young open star clusters) dominate 
the merger rate in the Universe, or whether multiple formation path- 
ways co-exist”. 

The majority of current NSC models yield a merger contribution 
of Rx O (10-1 Z 10!) Gpc™? yr! , and NSCs are observed to have 
atypical half-mass radius of ~ 0(10°) pc. As shown in Fig. 4, this would 
guarantee Naet ~ O(10) detections of SMBHBs after 7=10 years of 
observation. Inthe AGN channel, compact binaries merge even closer 
to the galactic centre, as they are within the extent d $ 0(10-') pc of 
the AGN disk. Thus, we find that current AGN models allow for 
Naet © O (107! — 10*)SMBHB detections. Note that our method allows 
SMBHB detections even if these channels only amount to a subdomi- 
nant fraction of the total rate of compact binary mergers’. Only in the 
least favourable scenario, in which all compact binaries are exclusively 
formed in the galactic field, do we expect no detections. In that 
case, we can assume the compact binaries to follow the galactic 
stellar mass distribution. Then, dis similar to the typical half-mass 
radius ~ 0(107) pc of massive galaxies, yielding Naet $ O(107!) 
for R $ 0 (10! — 10°) Gpc™° yr-!Hence, itis plausible that a substantial 
number of compact binaries experience a detectable modulation due 
to the GWs from central SMBHBs, unless all of them merge in the 
galactic field. 

Finally, our method can be used to distinguish different compact 
binary formation channels. For instance, the observation of a modu- 
lated compact binary would suggest a formation within a galactic 
nucleus, whereas the absence of such detection places strong 
upper limits of R $ O(10-3) Gpc * yr-KAGN) and $ 0(10-2) Gpc? yr 
(NSC) on the rate of compact binary mergers that originate from the 
galactic centres. 
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Fig. 4 | Total number and distribution of the expected SMBHB detections 
overaT = 10 years observation with a decihertz detector. Top: we show the 
expected number of SMBHB detections as a function of the compact binary 
merger rate R and the distance d from the Galactic centre at which the binaries 
typically merge. Compact binary mergers can occur through one or more 
formation channels. Rate predictions by several AGN*’ *°, NSC” and field 
channel models" ” are shown by coloured shaded areas and further detailed by 
the error bars in the top right. We assume a plausible range of 


ds 0(10-!) pc, ~ O(10°) pc and x O(103) pc to represent the typical half-mass 
radii of AGN discs, NSCs and large galaxies, respectively. Crucially, the large 
majority of models for the AGN and NSC formation channels would guarantee 
tens to hundreds of SMBHB detections. Bottom: we show the mass and redshift 
distributions of the detectable SMBHBs in terms of their probability density 
(PDF) and cumulative density functions (CDF). The method is sensitive to 
individual sources with total mass M = 10’-108 M, up to redshifts of z < 3.5. 


It is important to note that we have restricted our analysis to 
sources that exhibit at least one entire sinusoidal modulation. In this 
way, we exclude the possibility that the perturbation to the GW may be 
degenerate with any additional slowly accumulating phase shift. For 
instance, GWs from a compact binary orbiting a central (single or 
binary) SMBH may exhibit a detectable phase shift caused by Doppler 
motion if d < O(1) pc (ref. 39). This Doppler modulation may only be 
periodic, and thus degenerate with our effect, ifthe observational time 
is comparable or longer than the orbital period of the binary around 
the SMBH(B)s. For T= 10 years and M=108 M,, that is the case only if 
d < (10-2) pc. Such small separations arise exclusively in formation 
channels within migration traps of AGN disks*°“". Conversely, the simul- 
taneous detection of our effect and the acceleration of the compact 
binary” in the vicinity of a central SMBHBs would be an unambiguous 
signature of its existence. Furthermore, the periodicity of the sinusoi- 
dal modulation occurs on nanohertz frequencies. While strong gravity 
effects such as apsidal and spin precession may produce some oscilla- 
tory phase modulations, they occur on frequencies that are several 
orders of magnitude higher for binaries in the decihertz band. Addi- 
tionally, the latter is evolving throughout the chirp and can therefore 


not be degenerate with a constant low frequency modulation (see also 
Supplementary Fig. 1, where we tested the degeneracy explicitly with 
the effective spin parameter). 

To conclude, we have shown that the inherent sensitivity of pro- 
posed laser interferometers may be exploited to detect individual 
SMBHBsin the previously inaccessible mass range 210’ M,, at redshifts 
beyond the capacity of both current and future pulsar timing arrays. 
Currently, there are no other proposed methods that would realistically 
guarantee a direct detection of GWs from individual SMBHBs in the 
nanohertz band. A decihertz GW detector has exceptional scientific 
potential and should therefore be pursued with urgency, as it would 
open a large volume of the Universe to GW-based observations, simul- 
taneously in both the decihertz and the nanohertz band. 


Methods 

Parameter estimation 

A standard quantity to estimate the detectability of a GW strain h with 
a given detector is the SNR p defined as 


pih] = 4/ (hh), 6) 
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where we defined the noise-weighted inner product 


fox afb (f+ Of )a Sf) 
(ajb) -2f df a ; (4) 


Here, S,(f) is the noise power spectral density of the detector and the 
integration boundaries fin and fax correspond to the signal frequen- 
cies at the beginning and the end of the observation, respectively. In 
general, the waveform h depends on a set of physical parameters 
h=h(0°,6',@,...). 

Let ĝi be the correct values of the signal and 6! + A@ the best-fit 
parameters for a given realization of the noise. For a large SNR p, the 
measurement uncertainty AG follows a Gaussian distribution * 


1 


p(A6!) « exp (-5r)a000/), (5) 


where/;,is the Fisher information matrix defined as 


oh 
r=( 


agi (6) 


ze) 
dð) 


Thus, the root-mean-square error of 6' is 
(a0) = VV", (7) 


where J = [is the inverse matrix of T. Thus, the relative uncertainty of 
the ith parameter canbe estimated as V 5" /ĝ'. 

The limitations of the Fisher matrix formalisms have been thor- 
oughly explored in ref. 45. It is generally understood that, when it comes 
to GW parameter estimation, reliable estimates often require the full 
sampling of parameter posteriors to test for degeneracies and other 
features of the low-SNR regime. In GW parameter estimation, the goal 
is to maximize a likelihood £ of the form 


£ (O) œ exp [—(A(O) — Ain) |A(O) — AGinj))] 5 (8) 


where O; represent the parameters of an injected signal. In this work, 
we complement our Fisher matrix estimates with a series of 
MCMC-based parameter recovery tests. We then adjust the Fisher 
matrix detectability threshold to reflect the more conservative 
MCMC-based results. To explore the posterior distribution functions 
for the parameters O, we use the affine-invariant ensemble sampler 
emcee“ for MCMC proposed by ref. 47 and perform several tests with 
PN waveform models. In the MCMCs, we use uniform priors for all 
waveform parameters, except Amoa ANd finog for which we impose 
log-uniform priors between the wide ranges 107-10 * and10?-10* Hz, 
respectively. A realization of an MCMC test for a detectable modulation 
can be seen in Supplementary Fig. 1. In Supplementary Fig. 2, we show 
the sensitivity curve that results from injecting a grid of values of A moa 
and fmoa: This sensitivity curve is used in Fig. 2. 


Unmodulated waveform 

The simplest waveform from a circular binary can be obtained by 
modelling it as two Newtonian point masses m, and m, whose orbital 
frequency grows secularly due to the energy loss by GWs**: 


Fivewton(F) = SMSIf-76 expliv( FY] ©) 


wheref2 O, Dis the luminosity distance to the binary, Qis a numerical 
factor accounting for the angular emission pattern of the source and 
the antenna response of the detector to the GW, and there is a phase 


TU 


-5/3 
4 B 


U(f) (10) 


2nft, — Pe + 3 (sna 


Thus, this waveform can be written in terms of four physical 
parameters: an overall amplitude A = Qm*"/D, a chirp mass 
M = (mm)? (m, + m) , a collision time t, describing the point in 
time of the merger, and a constant phase @.. If the binary is at a cosmo- 
logical distance at redshift z, we need to replace its chirp mass by 
M > M(1+2). 

The advantage of Newtonian waveforms is that they may be easily 
treated in the Fisher matrix formalism and that the scalings and degen- 
eracies are well understood**. However, Newtonian waveforms do not 
suffice for the requirements of parameter estimation already in current 
GW data analysis pipelines. To surpass these limitations, we include 
PN corrections for our Fisher matrix estimate and the MCMC-based 
numerical tests. We include the full phase evolution for non-spinning 
binaries up to 3 PN, and include the effect of spin through an additional 
phase modification at 1.5 PN order. The coefficients of the GW phase are 
taken from ref. 48. We also include the effect of all higher-order modes 
up to powers of x75, where.x is the PN parameter, also taken from ref. 48. 
Effectively, this introduces also a dependency of the waveform onthe 
reduced mass 4 = m,m,/(m, + m,) and the effective spin parameter £. 


Modulated waveform 

We consider the situation where we have a ‘carrier’ source that emits 
GWs within the frequency bandwidth of our detector and another ‘back- 
ground’ source radiating at much lower frequencies. The background 
GW modulates the apparent frequency of the carrier GW by perturbing 
spacetime along the detector-carrier line of sight. Thus, the measured 
frequency can be written as” 


(11) 


1 
Simeasured(0) = TH + ATO x fO -SOKO 4 oO(@), 


where we substituted the unmodulated frequency f(t) =1/T(t) anda 
small modulation ¢(t) = AT(t)/T(t) < 1. Note that the modulating effect 
of the background GW applies to any generic periodic carrier signal. For 
instance, using telescopes to search for modulations of the rotational 
period of millisecond pulsars is the aim of large-scale observational 
campaigns to detect low-frequency GWs*°’°’. In this work, we assume 
that the carrier source is a binary emitting GW (see above) whose 
frequency evolves as’ 


—3/8 


= 1 fle -om 
KO = |k- om] (12) 
The modulation is given by” 
DiDi DiDi 
t= ——hA'T(t,x =O —— h(t- D,x = D), 13 
O 21 +b. D) y EXTD) 2(1+d. D) a ( xD Mes 


where x = Ois the position of the Earth, x = D the position of the carrier, 
x = bthe position of the background source, d = b - D the vector from 
the background source to the carrier, and h" the linear spacetime 
metric perturbation in the transverse-traceless gauge that is induced 
by the background source. Equation (13) is identical to the formalism 
used to describe the GW-induced modulation ofa pulsar®’, except that 
in that case Earth and pulsar are considered much closer than the 
SMBHB (D < d) so that the GW is incident from the same direction, that 
is, b x d. 

Inthis work, we are considering the opposite case where the carrier 
and the background source are close to each other, but far away from 
Earth, that is, D > dand b + d. In pulsar timing literature, the first term 
onthe right-hand side (r.h.s.) of equation (13) is usually referred to as 
Earth term and the second to as pulsar term. The amplitudes of the first 
and second term on the r.h.s. of equation (13) decrease with distance 
D and d, respectively. Due to the large distance to the background 
source, we can neglect the former in our work and focus on the 
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spacetime distortion the background source produces at the location 
of the carrier. Ifthe background source is an SMBHB witha chirp mass 
Mmoa that emits monochromatic GWs at a frequency fog, the distant 
modulation term can be written as*”” 


KO) = Amod cos(2Tfimnoat + Pmod); (14) 


where its amplitude is of the order Amod © (Whoa) Pusl3 id. Here, we 
averaged over a random orientation of d with respect to D and a random 
inclination of the orbital plane of the SMBHB. In the time domain, the 
GW strain from the carrier is given by* 


(384/5? 


2/3 t 
hÀ = p cos (on Smeasured(l’) ae) A (15) 


where M = m, + m, and u = m,m,/M are the total and reduced mass, 
respectively, and the orbital separation of the binary is 


256 1/4 


r(t) = (Flume) t-o (16) 


Considering the integral in equation (15), we have 
t t t 
WO = 2m f Freasea(t) de 2m f fede -2n f oae, a7) 


where the first integral on the r.h.s. evaluates to** 


T 5/8 


anf Roar = -2( 57) + Be, (18) 


where T=t,- tis the time to coalescence and the second integral 
becomes 


t 
Ahn fe, SAT 
an f zede ao o o o 
B,3 29, 2m2 (19) 
x [101 fnoaT sin (Pmoa) F, (S53. 2 16’ “J meat T ) 
+13 cos (Pmoa),F2 (=: > P ra + Pemod- 


The hypergeometric functions ‚F, evaluate to one in the limit 
(fmoaT ) > O where the period of the modulating GW is much longer 
than the observation time, and they are zero in the opposite regime 
(fmoaT ) > œ. Next, we consider the Fourier transform 


ACf) = J hoef de. (20) 
In the stationary-phase approximation, we get** 
Af) = 2 sif- -7/6 exp [i(2nft — Øf) — 1/4], (21) 
where @(f) = @(f(7)) and 
TF) = 5(8nf) P M-56. (22) 


Note that equation (21) reduces to equation (9) if there is no modulation 
Amod = 0. In general, the modulated waveform in equation (21) now 
depends on four additional physical parameters: an overall amplitude 
Amoa, a Modulating frequency fmoa and the phases @moq and Pe mod: 
However, we can further reduce the degrees of freedom to seven in 
total by omitting @.moq through an appropriate redefinition of 
Pc > Pe + Pcmoa- When we calculate the noise-weighted inner 
product in equation (4), the integration limits are chosen to be 


fin =f(T = 10 years) and fax = fisco = 1/12mV6(m, + m2) , that is, we 
assume that we observe all binaries during their last 10 years before 
they merge. In reality, when a detector begins to observe the binaries 
will be distributed uniformly in 7, that is, some of the binaries are in 
fact observed for less than 10 years before they merge 
(fmin > KT = 10 years)). Binaries that are observed for less than 10 years 
are less likely to exhibit one full cycle of a low-frequency modulation, 
meaning that the modulation is more liable to be degenerate with other 
environmental effects (see main text). However, for the typical carrier 
sources we are considering, more than ~90 % of the SNR is accumulated 
in the last year before merger. Indeed, we find that the bound on 
AAmoa/Amoa converges to the results shown in Fig. 2 well within a factor 
of 2 for sources that are observed for ~3 years or more. Thus, our 
population estimates will at most be reduced by ~30 % by this 
consideration. 


Carrier populations 
During the first three LIGO-Virgo-Kagra observing runs, the GWs of 
ninety mergers of BBHs, BNSs and neutron star—black hole (NSBH) 
binaries have been directly detected***®. From this sample, it is pos- 
sible to reconstruct the merger rate and parameter distribution of the 
underlying astrophysical population of compact binaries**. Owing to 
their stronger signal, most mergers have been detected of BBHs. Their 
inferred merger rate in the local Universe is*® 
Repu =17.9 — 44Gpc yr? (90% CN), (23) 
at a fiducial redshift z= 0.2. The rate is observed to scale with redshift 
as (1+z)*, where x = 2. on 7 for z<1. It is plausible to assume that this 
trend extends to redshifts z = 2-3 if it roughly follows the shape of the 
cosmic star formation rate density of massive stars”. In addition, the 
detection catalogue of the first three LIGO-Virgo-Kagra observing runs 
contains two BNSs and four NSBH binaries*®. The inferred merger rates 
in the local Universe are 
(24) 


Rens =10 — 1700Gpceyr-!_ (90% CN), 


Ryssen = 7.8 — 140 Gpc°yr™ (90% C1). (25) 


Due to the small number of detections these rate estimates are quite 
uncertain. Furthermore, no redshift evolution could be determinedin 
either cases. For the purpose of this work, we agnostically assume the 
same redshift scaling as for the BBH mergers. 

Regarding the chirp mass of the mergers, the BBHs follow a 
clumped distribution with pronounced overdensities around 
M =8.3+%} Mo and 27.9112 Mo. Meanwhile, the chirp masses of 
the six events involving a neutron star range from M = 11889005 Mo 
(GW170817 (ref. 68)) to 3.7 + 0.2 M, (GW190917 (ref. 69)). The chirp 
mass of the carriers affects our resultsi in two ways. On the one hand, a 
higher chirp mass increases the SNR of the carrier, which overall 
benefits the identification of a modulating background source. 
On the other hand, a higher chirp mass shortens the merger time and 
the time a binary spends in the frequency bandwidth of a detector 
(cf., equation (22)). This impedes the identification of a modulating 
background source in a more important way. To demonstrate this 
effect, we adopt M = 27.9 Mo from the first directly detected BBH 
merger (GW150914 (ref. 70)) as a fiducial value for the chirp mass of a 
rather massive equal-mass BBH and assume the chirp mass M = 1.188 Mo 
ofan equal-mass BNS as our standard choice for binaries with a neutron 
star component. 


Proposed decihertz detectors 

A GW detector that is sensitive to frequencies around = 0(107?)Hz is 
ideal for our type of analysis. In this frequency range, the compact 
binaries that we consider in this work spend atime T < f/f up to several 
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years before they merge in the bandwidth of LIGO-Virgo-Kagra. Moni- 
toring carrier binaries over several years allows us to detect frequencies 
of a modulating background source as low aS finog 2 1/T ~ O(107°) Hz, 
which would bea probe of the most massive black-hole binaries in the 
Universe. For this reason, we focus on existing detector proposals in 
the decihertz regime. For concreteness, we consider the space-based 
instruments DECIGO”’” and BBO” ” but stress that our analysis could 
be applied to any detector that is sensitive to those frequencies. We 
note that DECIGO is actively being developed, while we are not aware 
of any research efforts towards the realization of BBO. 

Both DECIGO and BBO consist of multiple (up to four) trian- 
gular constellations of spacecraft in heliocentric orbits, where the 
constellations are equally spaced around the ecliptic. DECIGO and 
BBO are designed to detect GWs at frequencies between 0.1 Hz and 
10 Hz such that they could provide sensitivity in a frequency range 
between the sensitive band of LISA and ground-based detectors such 
as LIGO-Virgo-Kagra. Both DECIGO and BBO are laser-interferometric 
detectors where the interferometer arms are rendered by beams of laser 
light that travel between the spacecraft. In general, the sensitivity of 
interferometers decreases above a frequency inversely proportional to 
the arm length, and therefore both DECIGO and BBO will employ arms 
shorter than those of LISA. DECIGO plans to use interferometer arm 
lengths -10° km; BBO’s arms are intended to havea length of -5 x 10* km 
(cf. LISA, which has a design with arm length ~2.5 x 10° km)?”*”°. While 
BBO has arms that are longer than those of DECIGO by a factor of ~10, 
DECIGO is designed with arms that contain Fabry-Perot cavities witha 
finesse of ~10, and the instruments therefore have almost equal strain 
sensitivities. The optical frequency response of both instruments is 
comparable as well: the sensitivity to GWs falls off linearly at frequen- 
cies above the pole frequency of the Fabry-Perot cavities in the arms of 
DECIGO, and similary above the light-crossing frequency of the arms 
of BBO. The design of DECIGO envisions the instrument to be limited 
by quantum noise across the entire band, specifically quantum radia- 
tion pressure noise and quantum shot noise. Extensive mitigation of 
classical low-frequency acceleration noise will be required to achieve 
this condition. Quantum radiation pressure noise falls off as f*, and 
quantum shot noise is constant at all frequencies; the crossover point 
of these two noise contributions is determined by the optical power 
circulating inthe arms. BBO is expected to be limited at low frequency 
by acceleration noise of the test masses due to a variety of sources, and 
at high frequencies by quantum shot noise. 

For both BBO and DECIGO, the triangular design of the instrument 
provides two independent interferometric measurements of incident 
GWs, which can be equivalently described as measurements from 
two independent L-shaped Michelson interferometers”*”””. Given the 
proposed experimental parameters, the quantum noise and optical 
response may be modelled, and the sensitivity of of each independent 
L-shaped interferometer can thus be described as an effective noise 
power spectral density S,. Here, we use the analytical expressions?’ 
for DECIGO 


2 
SPECIGO( F) = 7.05 x 10-48 jr (£) [r 


—4 
+4.8 x 10-5(-£) 1 yz? (26) 
1Hz (£) 
fp 
—4 
+5.33 x 10-?(-) Hz, 
1Hz 
where f, = 7.35 Hz and for BBO 
2 
BBO = -49( f -1 
$880 f) = 2.00 x 10 (4) Hz 
+4.58 x 10-49 Hz (27) 


+1.26 x 10-3(£) “Ha 


The sensitivity further depends on the implementation of the instru- 
ments; the achievable SNR increases as Niro, Where Mo is the number 
of independent L-shaped Michelson interferometers that the detector 
furnishes. In this paper, for both DECIGO and BBO, we take a sensitivity 
corresponding to one triangular constellation, such that N fo = 2, that 
is, we divide equations (26) and (27) by V2. Furthermore, we assume a 
nominal mission lifetime for both instruments of T=10 years. In assess- 
ing the sensitvity of DECIGO/BBO, we also consider the angular 
response factor Q (cf. equation (9)). Averaging over the antenna pattern 
for sky positions (9, @) ofan L-shaped Michelson and all possible binary 
orbital inclinations (1) gives” 


Via. @.0P) = F, 


which we will use for Qin this work. 


(28) 


Background source population 

The background sources considered in this work consist in individual 
SMBHBs, which can span over several decades in mass, mass ratio and 
orbital separation. The SMBHB merger rate is highly uncertain, and 
constraining it is one of the major scientific goals of proposed GW 
observatories such as LISA®”, TianQin”’, pulsar timing arrays” ° and 
other electromagnetic searches****. For the purposes of this work, we 
broadly follow the procedure detailed in ref. 86 and link the SMBHB 
merger rate to the well-established halo merger rate, based on the two 
Millenium simulations”: 


Mhato a 2 ç a b4 
ots) &: ovo |(5) Jara ; 


where His the total number of mergers that a halo of mass Mpa experi- 
ences over cosmic time, €< 1 is the halo merger mass ratio and 
the best-fit parameters are given by [B,, B, b,, b», b>, b4] = 0.0104, 9.7 
2 x 9.72, 0.133, - 1.995, 0.263, 0.0993]. The SMBHB merger rate W,,can 
be obtained by multiplying the halo merger rate with the SMBH mass 
function 


2 
d*H 6,( 


oH aB (29) 


aN = P(Mpal team?) OE 
= alo» 


aM.dédz dl + 2 dM. dédz 630) 


(€ Zaeri), 


where c = 1and where we must supply an occupation fraction P, a rela- 
tion between SMBH and halo masses, and a delay prescription between 
the nominal halo merger and the actual SMBHB merger. Considering 
that the most easily detectable modulations will be naturally produced 
by heavy, low-redshift binaries we can tailor equation (30) to SMBH 
masses of M. > 108 M, within moderately low redshifts z < 4. First, we 
adopt the SMBH mass function reported in ref. 88, which we denote as 


dn. 


———_., 31 
dlog,,M. 6) 


ġ= 


where @ is the symbol customarily used to denote the SMBH mass func- 
tion in units of Mpc” yr“ dex”. Secondly, we adopt a simple relation 
between halo and SMBH mass from ref. 89: 


3/2 


Mpato(1 +z) Mo. (32) 


2x 107 Mo 


which is broadly consistent with both theoretical and observational 
constraints”®”". Given the selection effects towards higher masses, 
we can safely assume an occupation fraction P = 1 as it is strongly sug- 
gested by simulations” and observations of AGN”. Finally, the time 
delays between the halo merger and the eventual SMBHB merger 
are expected to be of the order 108-10° years for the SMBHBs we are 
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considering’*”. As seen in Supplementary Fig. 2, we have tested the 
effect of introducing constant time delays of up to 1 Gyr by shifting the 
merger redshift to z > z- Az(t), where Tis the time delay and we assume 
astandard cosmology. We find that our results do not change appreci- 
ably, decreasing by at most ~20% and only at redshifts higher than ~3. 
This result runs counter to the common expectation that introducing 
delays reduces the number of detected massive BH mergers for GW 
detectors’®”. In our work, it can be attributed to two factors. Firstly, 
the halo merger rate reported in the Millennium simulation only weakly 
depends on redshift, scaling as ~(1 + z)”. Secondly, in our setup, shifting 
SMBHB mergers to lower redshift automatically increases the SNR of 
the detected stellar-mass compact object binaries while simultane- 
ously reducing their cosmological redshift. This has the consequenceto 
greatly expand the range of detectable modulations in both frequency 
and amplitude, counteracting the diminished number of SMBHB 
mergers. The fact that these effects seem to cancel out each other is 
somewhat of a coincidence, which we attribute to the explicit form of 
equation (29). 


Rate estimate 

Given a population of both carriers and background sources, we can 
estimate the mean number of compact binaries that could exhibit a 
detectable modulation. To achieve this, we must select all SMBHBs that 
are emitting loud GWs in the appropriate frequency range for the 
sensitivity curves detailed above (cf. Fig. 2). The first step is to distribute 
the SMBHB over several frequency bins, by replacing the time derivative 
with an integration over fnod: 


N.. fnat 


Note that, since we are considering SMBHBs in the GW-driven regime, 
we can simply use equation (12) to describe the frequency evolution. 
Then, we can integrate the differential contributions over frequency 
and mass ratio, using a Heaviside function Oto only select SMBHB that 
would produce detectable modulations, as defined by acertain confi- 
dence threshold dget: 


d2Nmod d*N AA 
a EN on ra) n mes ) d dé. 
dM, dz / f dM. dédfmnod dz (ae Amod Vinod £ 


The Heaviside function depends on the sensitivity curves detailed 
above (cf. Fig. 2): 


dN.. ; 


dN. 
dt = 
df, roa ™8 


dfmod 


(33) 


dfmoa 


(34) 


O = OZ, Amod mod: M, tes Pes Omoa)» (35) 


where we initially choose a threshold of AA mod = Amog that is, dget = 1, 
to represent a SMBHB modulation, the parameters of which can be well 
constrained. Note that, typically, constraining the parameters of a GW 
to lodoes not constitute a sufficient threshold to avoid the possibility 
of false signals. In our case however, we are considering modulations 
that affect GWs signals with SNR of order few 107, and it is not immedi- 
ately clear whether the same intuition applies, nor if the posterior 
distributions of the parameters are Gaussian. To test this, we have run 
additional MCMC tests in which the null hypothesis, that is, the absence 
ofa signal, is tested against several realizations of Gaussian noise. We 
observe that the MCMC walkers freely sample regions in which the trial 
amplitude is very low. However, they are strongly discouraged to sam- 
ple large values of Amoa, Which are inconsistent with the absence of a 
signal. The posterior probability for the recovered modulation ampli- 
tude Amoq Sharply decreases above a certain threshold Ayp associated 
to the null hypothesis 


for Amod < Anu(Smoa) 


P(Amod) ~ 
mo for Amod > Annul fmoa). 


const. 
l (36) 


More precisely, the boundary is a half-Gaussian, for which the inflec- 
tion point defines the loy,, confidence threshold between the 
allowed and disallowed regions. We find that the boundary at which 
the null hypothesis can be confidently excluded always lies about 
one order of magnitude below our choice of sensitivity curve, espe- 
cially at lower modulation frequencies, where the majority of detec- 
tions take place (Supplementary Fig. 2). In this sense, our choice of 
using alo boundary simply refers to the desired accuracy in the 
recovered modulation parameters, though it is likely that a more 
complex noise model would affect this conclusion. Therefore, the 
lo detection criterion of Amoa only defines how accurate we could 
determine its value when a modulation is present and does not imply 
a large false positive probability at which we would wrongfully 
conclude the existence of a modulation. We have additionally tested 
how the number of detection scales as a function of the desired 
posterior widths, finding that 


Naet(Odet) & oz. (37) 


This means that a substantial fraction of the detections will result in 
better-constrained posteriors. 

The next step is to distribute the total amount of compact object 
mergers into their host galaxies and across cosmic time. To enable 
this calculation, we require to know the differential compact object 
binary rate per halo. The most natural assumption is that the latter 
must be proportional to the stellar mass of the galaxy that occupies 
ahalo of a given size. Interestingly, for the high-mass galaxies likely 
to host heavy SMBHs, the stellar mass—-halo mass relation has a turno- 
ver point at a value of ~1/100, varying by approximately one order of 
magnitude over the large range of 10" M, < Mhaio < 10% Mo. For the 
purposes of this work, we use a fit of the stellar mass to halo mass 
ratio F* (refs. 98,99): 


-1 


F* (Mhraio) = 2D,(1 + 2)" (“ee ) (“i )] (8) 
2 2 

Dy = 108+? Mo, (39) 

B = 26, + 6s, (40) 

y=6,(1+2)”, (41) 


where the best-fit parameters are (D,, 6,, ô», 63, 64, 65, 6, 67) = (0.046, 
-0.38, 11.79, 0.20, 0.043, 0.96, 0.709, —0.18). Thus, the number of 
compact object binaries is distributed as 


dR R dha 
— <M = EOM F* (Mato) Malo —— + 


= (42) 
dM, halo N dM, halo 


where the normalization is given by an integral over the total stellar 
mass contained in all halos in the Universe: 


dn 
N= J = Maora See Mt (43) 


dMhalo 


and we obtain the halo mass function by inverting equation (32) 
and adjusting the differentials accordingly. Note that relaxing the 
assumption of a constant occupation fraction does not greatly affect 
the normalization because a large fraction of the total mass budget is 
composed by massive halos above -10° M,. 

To estimate the probability of a compact object merger taking 
place in a galaxy that simultaneously hosts a modulating SMBHB, we 
suppress the former’s rate by a further factor representing the number 
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of halos that contain such an SMBHB with respect to the total number 
of halos, ina given redshift shell with volume dV,: 


(44) 


PREM _ Room d2nmod ANhato dV, = 
AMhato dz ~ dMnaiodz ae aie =) , 
where once again we make use of the halo mass-SMBH mass relation 
and adjust the differentials accordingly. Note that here the compact 
object merger rate must be transformed from the local source frame 
to the observer frame, by dividing the observed redshift dependence 
by a factor (1+ z) and expliciting the differential d/dz. 

Finally, the total number of compact object mergers with a detect- 
able modulation N™°4 can be obtained by integrating and multiplying 


COM 
by the total observation time: 


d2.pmod 
AMpato dz = Nadev 


Nmod =T 
AMpalo dz 


COM (45) 


where, equivalently, Ngo, isthe number of indirectly detected SMBHB. 
Note that in reality the number NT is actually the expectation 
value of a process that involves sampling both compact object binaries 
and SMBHB from a large pool of halos, analogous to a binomial distribu- 
tion, excluding the unlikely possibility of having two separate compact 
object binaries merging simultaneously with an SMBHB in the same 
galaxy. Therefore, it carries an intrinsic uncertainty approximately 


proportional to its square root. 


Data availability 
The data and codes for reproducing the results of this work are avail- 
able via GitHub at https://github.com/stegmaja/GW-GW-Modulations. 
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